#include <fstream>
#include <vector>
#include <string>
#include <iostream>

#include "NURBS_Bezier.hpp"
#include "NURBS_Tools.hpp"
#include "NURBS_FileIO.hpp"
#include "DArray.hpp"
#include "DMatrix.hpp"
#include "BernsteinBasis.hpp"
#include "MeshProcess.hpp"
#include "MeshInfoProcess2D.hpp"
#include "QuadElement.hpp"

using namespace std;

int main(int argc, char *argv[])
{
  /* 
     vector<double> knot(1003, 0);
     knot.push_back(0.0);
     knot.push_back(0.0);
     double h = 0.001;
     for(int i=2; i<1001; i++)
     knot.push_back((i-1)*h);
     knot.push_back(1.0);
     knot.push_back(1.0);



     int cDegree = 3.0;
     int dim = 3;
     vector<NURBS_T::BezierElem*> bSeg;
     vector<double> ctrlPts;
     vector<double> bCtrlPts;
     NURBS_T::decomposeCurve(knot, cDegree,ctrlPts, dim, bCtrlPts, bSeg, false );

     for(int i=0; i<4;i++)
     {
     for(int j=0; j<4; j++)
     cout<<bSeg[1]->coefficient(i,j)<<'\t';
     cout<<endl;
     }

     for(int i=0; i<(*bSeg[2]).gPtsMap.size();i++)
     cout<<(*bSeg[2]).gPtsMap[i]<<endl;

     vector<NURBS_T::BezierElem*>::iterator it;
     it = bSeg.begin()+1;
     for(it=bSeg.begin(); it<bSeg.end();it++) 
     (*it)->~BezierElem();

*/

  /*
     DArray test1(10), test2(10);
     test1[1] = 1.0;
     test1.print();
     test1*=2.5;
     test1.print();
     test2[1] = 1.2;
     test2[4] = 1.9;
     test1+=test2;
     test1.print();
     */




  vector<double> sKnots, tKnots, uKnots;
  int sDegree, tDegree, uDegree, numCPts;
  vector<double> ctrlPts;
  int nqpx = 3;
  int nqpy = 3;

  ifstream infile( "geometry.txt", ifstream::in );

  NURBS_T::readFile( infile, sKnots, tKnots,uKnots,
      sDegree, tDegree, uDegree, numCPts, 
      ctrlPts);

  NURBS_T::projectUp( ctrlPts, 3 ); 

  vector<double> insertKnots1, insertKnots2;
  double h = 0.2;

  insertKnots1.push_back(0.2);
  insertKnots1.push_back(0.3);
  insertKnots1.push_back(0.4);
  insertKnots1.push_back(0.6);
  insertKnots1.push_back(0.8);

  insertKnots2.push_back(0.2);
  insertKnots2.push_back(0.4);
  insertKnots2.push_back(0.5);
  insertKnots2.push_back(0.6);



  NURBS_T::knotRefinementSurface( sKnots, sDegree, 
      tKnots, tDegree,
      4,ctrlPts, insertKnots1, 's' );

  NURBS_T::knotRefinementSurface( sKnots, sDegree, 
      tKnots, tDegree,
      4,ctrlPts, insertKnots2, 't' );


  vector<NURBS_T::BezierElem*> seg_x, seg_y;
  vector<double> fake_ctrl;
  vector<double> fake_bCtrlPts;
  NURBS_T::decomposeCurve(sKnots, sDegree, fake_ctrl, 4, fake_bCtrlPts, seg_x, false );
  NURBS_T::decomposeCurve(tKnots, tDegree, fake_ctrl, 4, fake_bCtrlPts, seg_y, false );


  

  MeshInfo_2d mesh2d(sKnots, tKnots, sDegree, tDegree, nqpx, nqpy);

  

  vector<double>::iterator ita;
  for(ita = sKnots.begin(); ita<sKnots.end();ita++)
    cout<<*ita<<'\t';
  cout<<endl;

  for(ita = tKnots.begin(); ita<tKnots.end();ita++)
    cout<<*ita<<'\t';
  cout<<endl;




  QuadElement test( 2, mesh2d, ctrlPts, 4, seg_x, seg_y);  
  test.ElemShapeRoutine(0.5, 0.5);

  for(int i=0; i<9; i++)
    cout<<test.R[i]<<endl;
  cout<<endl;






  /* === Clean the dynamic array === */ 
  /* 1. Extraction operator */
  vector<NURBS_T::BezierElem*>::iterator it;
  for( it = seg_x.begin(); it<seg_x.end(); it++ )
    delete *it;
  
  for( it = seg_y.begin(); it<seg_y.end(); it++ )
    delete *it;




  /*
     NURBS_T::knotRefinementVolume( sKnots, sDegree, 
     tKnots, tDegree,
     uKnots, uDegree,
     ctrlPts, insertKnots, 'u', 4 );
     */

  /*

     vector<double> outDisp;

     int ii;
     for( ii=0; ii<ctrlPts.size()/4; ii++ )
     {
     outDisp.push_back(0.0);
     outDisp.push_back(0.0);
     outDisp.push_back(0.0);
     outDisp.push_back( ctrlPts[ii*4+3]);
     }

     outDisp[8] = 1.0;
  //outDisp[0] = 1.0;
  outDisp[4] = 1.0;
  NURBS_T::projectUp(outDisp, 3);

  cout<<"The control points are:"<<endl;
  for(int i=0; i<ctrlPts.size()/4;i++)
  {
  cout<<ctrlPts[i*4]<<'\t'<<ctrlPts[i*4+1]<<'\t'<<ctrlPts[i*4+2]<<'\t'<<ctrlPts[i*4+3]<<endl;
  }

  cout<<"The output disp are: "<<endl; 
  for(int i=0; i<ctrlPts.size()/4;i++)
  {
  cout<<outDisp[i*4]<<'\t'<<outDisp[i*4+1]<<'\t'<<outDisp[i*4+2]<<'\t'<<outDisp[i*4+3]<<endl;
  }
  */

  /*
     NURBS_T::writeSurface( "resultSurface", sKnots, sDegree,
     tKnots, tDegree, ctrlPts, outDisp, 3, 100 );


     NURBS_T::projectDown(ctrlPts, 3);


     vector<double> dummy;
     dummy.push_back(0);



     NURBS_T::writeControlNet( "outNet", sKnots, tKnots, dummy, sDegree, 
     tDegree, 0, ctrlPts, 3 );
     */


  return 0;
}
